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Abstract 

Moments of the quark density, helicity, and transversity distributions are cal- 
culated in unquenched lattice QCD. Calculations of proton matrix elements 
of operators corresponding to these moments through the operator product 
expansion have been performed on 16 3 x 32 lattices for Wilson fermions at 
(3 = 5.6 using configurations from the SESAM collaboration and at f3 = 5.5 
using configurations from SCRI. Onedoop perturbative renormalization cor- 
rections are included. At quark masses accessible in present calculations, 
there is no statistically significant difference between quenched and full QCD 
results, indicating that the contributions of quark-antiquark excitations from 
the Dirac Sea are small. Close agreement between calculations with cooled 
configurations containing essentially only instantons and the full gluon con- 
figurations indicates that quark zero modes associated with instantons play 
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a dominant role. Naive linear extrapolation of the full QCD calculation to 
the physical pion mass yields results inconsistent with experiment. Extrapo- 
lation to the chiral limit including the physics of the pion cloud can resolve 
this discrepancy and the requirements for a definitive chiral extrapolation are 
described. 
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I. INTRODUCTION 



The quest to understand the matter of which our universe is composed will remain funda- 
mentally incomplete until we understand how the quark and gluon structure of the nucleon 
arises from QCD. The nucleon has many remarkable properties. Because of confinement, the 
quark-gluon structure of hadrons differs essentially from that of any other known composite 
systems. Gluons in QCD are essential dynamical degrees of freedom, unlike the boson fields 
in atoms and nuclei which may be subsumed into a two-body interaction thereby reducing 
these systems to purely fermionic degrees of freedom. Almost all of the mass and approx- 
imately half of the momentum and angular momentum of a nucleon arises from gluons. 
Indeed, even the net spin 1/2 arises from a rich and complicated combination of orbital 
and intrinsic angular momentum of the quark and gluon fields. Since the usual analytic 
tools of theoretical physics have proven inadequate to solve nonperturbative QCD, the only 
known way to solve, rather than model, QCD is numerical solution of lattice field theory. 
Hence, the ultimate goal of this work is to use lattice QCD to understand the structure of 
the nucleon. This understanding has two important but distinct aspects. 

One aspect of using lattice QCD to understand nucleon structure is the quantitative ab 
initio calculation of experimental observables. Since the experimental discovery of quarks in 
the nucleon over a quarter of a century ago, there has been a huge investment internationally 
in using high energy scattering to measure the light cone distribution of quarks and gluons 
in the nucleon. As a result of several decades of experimental effort at SLAC, Fermilab, 
CERN, and DESY, we now have a detailed knowledge of the quark density and helicity 
distributions and of the gluon distribution £[]-§]. in addition, major new experiments are 
being planned at these facilities as well as at Jefferson Lab and RHIC to map out the quark 
and gluon structure of the nucleon in even more detail. Now that the techniques of lattice 
field theory and computer technology have developed to the point that it will be possible 
to solve QCD with a precision comparable to experimental measurements, it is essential to 
complement this massive experimental investment with a commensurate theoretical effort in 
lattice QCD. Because deep inelastic lepton scattering measures correlation functions close 
to the light cone, structure functions are intrinsically Minkowski and cannot be calculated 
directly in lattice QCD. However, using the operator product expansion, it is possible to 
calculate their moments. In this work, we report the first calculations of these moments 
in full QCD ||. As will be emphasized below, these present calculations are subject to 
significant limitations due to current computer resources. 

The second aspect is understanding the basic mechanisms underlying nucleon structure 
- that is, how QCD actually works. Hence, we seek to use the lattice as a tool for insight 
as well as for numbers. There are several ways lattice calculations can provide insight. 
One is to calculate the overlap between a trial wave function and the exact nucleon wave 
function to explore the role of various degrees of freedom and variational parameters. Thus, 
this work will present a prototype variational calculation in which the rms radius of a trial 
function is varied. A second is to study the contributions of different classes of Feynman 
diagrams that correspond to separate lattice contributions to operators. An example is the 
class of connected diagrams considered in this work. Finally, since the lattice Monte Carlo 
calculation stochastically samples gluons distributed according to the QCD action, one can 
identify and study the dominant configurations. Hence, this work will study the role of 



3 



instantons and their associated zero modes in calculating moments of structure functions. 

There have been several calculations of moments of structure functions in recent years, 
including contributions of connected diagrams to low moments of the spin independent and 
longitudinal spin dependent structure functions in quenched lattice QCD [|10| p72|,|Qf| , discon- 



nected contributions for the axial and tensor charge and the axial charge in full as well 
as quenched QCD [fLlflj . Interestingly, in contrast to spectroscopy, there are significant dis- 
crepancies between these hadron structure calculations and experiment. Whereas quenched 
masses for light hadrons are typically accurate at the 5% level, the axial charge is typically 
10 - 20 % low ]T5| and the first moment of the spin averaged structure function is of the 
order of 50% high |10|,|TTJ. By the variational principle, we know that an error of order e in 
the wave function only produces an error of order e 2 in the expectation value of the energy, 
or mass in this case, so it is consistent that these other observables should be much better 
diagnostics than masses of the errors in the lattice calculation of hadron structure. Hence, 
we need to identify and correct the source of these substantial discrepancies. 

There are three major limitations in structure calculations imposed by limitations in 
computational resources. The first limitation is errors in approximating the continuum limit. 
The approach to the continuum limit has been studied systematically in quenched QCD by 
the QCDSF collaboration |Tt]] using Wilson and clover improved actions and extrapolating 
in inverse coupling (3. A particularly careful study was carried out for the axial charge which 
showed that even when the f3 dependence is included, there remains a 10% discrepancy with 



experiment 15 



The second limitation is the quenched approximation, which ignores the contribution of 
dynamical quark-antiquark excitations in the Dirac sea. It is now known that whereas full 
QCD gives a good description of masses of hadrons containing strange quarks, quenching 
causes significant discrepancies in masses. Physically, this is consistent with the fact that 
in quenched QCD the coupling runs too fast, producing a potential that becomes too weak 
at the short distances probed by the relatively heavy strange quark. Although there is no 
simple argument suggesting the sign of the resulting discrepancy in other specific observ- 
ables, it is important to see to what extent quenching causes the discrepancy in the axial 
charge and the first moment of structure functions. Thus, a primary goal of this work is to 
study the role of the quenched approximation. Although the effects of the lattice spacing 
and lattice volume are still significant potential sources of error, it is meaningful to compare 
quenched and unquenched calculations on comparable lattices in physical units to explore 
the magnitude of errors due to quenching. Hence, we will compare full QCD calculations 
using configurations produced by the SESAM collaboration [16[ at /? = 5.6 with quenched 
calculations on comparable lattices. Whereas it has been conjectured that the discrepancy 
between lattice calculations and experimental moments of structure functions arises pri- 
marily from quenching, one major result of this work is to show that at the quark masses 
attainable at present, this explanation is wrong. 

The final limitation is errors in calculating the large volume, small quark mass limit. 
Since the lattice volume must be large enough that the pion Compton wavelength fits well 
inside and inversion algorithms become less efficient for light quarks, it is presently impossi- 
ble to perform full QCD calculations with dynamical quark masses corresponding to physical 
pions or even to pions sufficiently light that chiral perturbation theory gives reliable extrap- 
olations. Thus, a major uncertainty is the extrapolation from the relatively heavy quark 
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masses used in lattice calculations to the small quark mass required to produce the physical 
pion. This is particularly worrisome for calculating hadron structure because of the major 
role played by the pion cloud in the nucleon. Since a small box and heavy quarks suppress 
the pion cloud, it is not surprising that the nucleon magnetic moment, much of which comes 
from the pion current, is low or that the axial charge is too small. Indeed, model estimates 



17|Jlql show that both discrepancies plausibly arise from omission of the full contribution 
of the pion cloud in present lattice calculations. Hence, a second major result of this work 
is to show how an extrapolation incorporating the leading effects of chiral symmetry can 
simultaneously resolve the discrepancy in the lowest three moments of the spin averaged 
structure function. 

The outline of this paper is as follows. Section II provides the background and defines 
the operators we evaluate. Perturbative renormalization is discussed in Section III and tech- 
nical details of the lattice calculation are given in Section IV. Section V presents the results, 
including the overlap of trial wave functions with lattice hadron ground states, comparison 
of observables calculated with quenched and unquenched configurations to study quenching 
errors, comparison of observables calculated with cooled and uncooled unquenched configu- 
rations to study the role of instantons and their associated zero modes, and comparison of 
linear and chiral extrapolations with phenomenology. A summary and conclusions are given 
in Section VI. 



II. BACKGROUND 

A. Moments of nucleon light cone quark distributions 

By the operator product expansion, moments of the linear combinations of quark and 
antiquark distributions in the proton 

( x ») q = f 1 dxx n (q(x) + {-l) n+1 q(x)) (1) 



x n ) Aq = t dxx n (Aq(x) + (-l) n Aq(x)) 
Jo 



(x n ) 5q = / dxx n (5q(x) + (-l) n+1 5q(x)), 
Jo 

where the quark density, helicity, and transversity [Bl distributions 



Q = <?t + Qi (2) 
Ag = g T - qi 

5q = q T - q ± , 

are related to the following matrix elements of twist-2 operators 

2 (x^X P,* ■ ■ ■ P, n = 2 Y,(PS\ (j) ri^D,, ■ ■ ■ D, n} V\PS) (3) 

2 T (x n )A qr s {a p m ■ ■ ■ p, n} = -{ps\ Q) • ■■D, n} r\ps) 



n + 

2 (x n ) 5qr ^P W P W • • • P, n} = (PS\ fi ) pw^Dn ■ ■ ■ D, n} V\PS) 



m N 
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Here, D = D — D, r denotes the quark flavor, x denotes the momentum fraction carried by 
the quark, S 2 = m 2 N , { } and [] denote symmetrization and antisymmetrization respectively, 
and the mixed symmetry [ { ] } term is first symmetrized and then antisymmetrized so that 
it is written explicitly as 



O 5 

[c{Ml]-Mn} _ n _|_ 



- (n 5 _ n 5 _l n5 — O 5 -\ ^ (4) 



We note that the odd moments (x n ) q are obtained from the spin-independent structure 
functions Fi or F 2 measured in deep inelastic electron or muon scattering 

j o 1 dxx n - 1 F 1 (x,Q 2 ) = \c°JC?ltf) e 2 r {x n ~ l )M (5) 
f dxx n ~ 2 F 2 {x,Q 2 ) = C^Q'/fi 2 ) j: e 2 r {x n ' l )M, 

JO r 

and even moments of {x n )& q are determined from the spin-dependent structure function gi 
f Q dxx n g x {x,Q 2 ) = \c a n (Q 2 /n 2 ) Y e 2 r 2(x n ) A M, (6) 
where e r is the quark's electric charge, and C n denotes the Wilson coefficient. Note that the 



moments (x n ) q and {x n )& q are proportional to the quantities v n+ \ and a n defined in Ref. JlC 

{x n ) q = v% (7) 

In addition, the two spin-dependent structure functions g\ and g 2 also determine the 
quantity d n 

— d r n S [a P M ■ ■ ■ P, n} = -{PS\ (-) ri 5 JwD M ■ ■ ■ D, n] V\PS) (8) 
which is a twist-three operator and does not have a simple interpretation in terms of parton 



distribution functions pTJ. However, since with Wilson fermions, 757(0- • • ■ D^} mixes 

with the lower dimension operator ^757[o- 7{ w ] • • • D/j, n }, it is not possible to compare with 
phenomenological results using the perturbative renormalization constants and mixing co- 
efficients calculated in this work. Either nonperturbative renormalization is required as 
has been carried out in []55|], or the operators need to be recalculated with overlap fermions 
P3[ or some alternative formulation for which mixing with lower dimension operators does 
not occur . 

Even moments {x n ) q are obtained from deep inelastic neutrino scattering, and in addition, 
a variety of other processes have contributed to what is now a detailed empirical knowledge of 
the quark and antiquark distributions in the nucleon. Hence, we will subsequently compare 
our results with moments calculated from the CTEQ, GRV, MRS, GRSV, and GS and 
global fits to the world supply of data |j]-|5|]. The moments of parton distributions (x n ) (/z), 
(x n )Aq(n), and (x n )s q (fi) are scheme and scale dependent, and we will convert our lattice 
matrix elements to the MS scheme and evaluate them at the scale /1 2 = -y ~ 4 GeV 2 . 
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B. Lattice operators 



Our objective is to calculate matrix elements of traceless, appropriately symmetrized 
and antisymmetrized operators of the general form 

= VSr 7M1 B M2 ---B Mn v, (9) 

where T = 1,75, or 7570-, on a hypercubic lattice to approximate the corresponding contin- 
uum operators as accurately as possible. Hence we choose representations of the hypercubic 
group, H(4) to eliminate operator mixing as much as possible, and after fulfilling this 
objective, to minimize statistical errors by including as few nonzero components of the nu- 
cleon momentum as possible. 

Since H(4) is a subgroup of the Lorentz group, irreducible representations of the Lorentz 
group are in general reducible under the H(4) group, and we choose the representation to 
optimize the approximation. It is essential to choose a representation that does not mix with 
lower dimension operators, since the coefficients would increase as l/a n in the continuum 
limit. In addition, because of the possible inaccuracy of perturbative mixing coefficients 
and the difficulty of determining mixing coefficients nonperturbatively, it is desirable to 
avoid mixing with operators of the same dimension as well. In choosing between operators 
with the same mixing properties, it is desirable to use a nucleon source with as few non- 
zero spatial momentum components as possible, since each projection introduces substantial 
stochastic noise. In fact, we will subsequently show that with available configurations, it is 
not possible to obtain adequate statistics in any momentum sector other than p = 0. Since 
any expectation value of an operator with tensor index j is proportional to Pj (or Sj for spin 
dependent quark distributions), the nucleon must have an additional momentum component 
projection for each new distinct tensor index that is added to an operator. Hence, the goal 
is to limit the number of distinct spatial indices. Eventually, as one proceeds to higher 
moments of quark distributions, all the space-time indices are exhausted and it becomes 
impossible to avoid mixing with lower dimension operators. 

The representations we have chosen for our operators using these criteria are enumerated 
in Table |I[ To illustrate the selection process, we describe selection of the spin-independent 
operators and analogous analysis yields the remaining operators. 

To measure (x) q one needs to calculate matrix elements of the traceless part of the 

operator q r ){ ll Dv}q ) which belongs to the representation (1,1) in the continuum decomposition 

® = (o, 0) © (1, 0) © (0, 1) © (1, 1). 

On the lattice, the nine dimensional representation (1,1) splits into two irreducible repre- 
sentations, 3^ and 63", both of which are symmetric and traceless, where the notation for 
representations is described in the caption of Table |H[ As a consistency check, it is desirable 
to calculate operators from each representation. For the first operator, denoted (x)^\ we 
select the basis vector of 63": 

(P\W{iD4 } q\P) = 2(x)V>-P {1 Pi h 
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and for (x)^>, we choose the basis vector of 3^: 

(P|g7 4 B 4 g - ^(qiiD iq + qj 2 D 2 q + qi£*q)\P) = 2(x)® ■ (P 4 P 4 - \p 2 )- 

Note that since (x)^ involves 74D4, it can be measured with P = whereas since (x)^ 

involves 71 D4, it requires a state projected onto non-zero Pi. 

For (x 2 ) q , none of the three (symmetric) representations 4j~ is appropriate, since they 
are not traceless and hence mix with lower- dimensional operators. The only representations 
with two distinct indices are the one 8^, which is not symmetric and must therefore be 
rejected, and the two 8^f 's 

q (1W1D1 - ^(74P>2P>2 + 74P>3P>3)) q, 

and 

q (71B4B1 + 71B1B4 - - (12DW2 + 72B2B4 + 73P>4P>3 + rsDsD^ q. 



which mix as discussed in Refs. [p5| -p7| . 

For {x 3 ) q the following representations have positive charge conjugation and do not mix 
with lower dimensional operators: 2f, 2%, 3^, 3%, 3%, 62" and 64 . However, the only 
representations that require a single non-zero momentum component are the two which 
generically could mix with each other but do not mix at the one loop level for Wilson or 
overlap fermions PBL|2"3L 



Note that in addition to the mixing discussed above, in full QCD there is also mixing 
between gluonic operators and flavor- singlet fermion operators for moments of the quark 
density and helicity that will not be considered in this work because we have not yet evaluated 
lattice matrix elements of the relevant gluon operators. 



III. PERTURB ATIVE RENORMALIZATION 

Since the phenomenological light cone quark distributions with which we compare our 
lattice results are extracted from experimental data using the MS renormalization scheme, 
we have converted our lattice calculations to the MS scheme in 1-loop perturbation theory 
using 

0F(Q 2 ) = £ fa + -A ^ (lW MQV) - {Bff* - Bp))) ■ Of ATT (a 2 ). (10) 



j 



The anomalous dimensions 7^ and the finite constants B^ are given in Table |TXT| for Wilson 
fermions and the specific operators considered. 

The Z factors that convert lattice results to the MS scheme at scale Q 2 = 1/a 2 are equal 

to 

Z(g 2 = 6/0) = 1 - ^ \ (B^ - , (11) 
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and are tabulated for two typical values of (3. Note that all of the moments of quark 
distributions calculated in this work will be presented at the scale of fi 2 = 4 GeV 2 in the 
MS scheme. 

Details of perturbative renormalization may by found in Refs. |29|j30|j2^| . The results 



in Table |TTT| are taken from Refs. in which the renormalization factors of (x)$ q and 

di for Wilson fermions were calculated for the first time and the remaining operators were 
checked with earlier results in Refs. J5T] , p3| JT0| , |2B| J5^] , revealing a discrepancy in the case 
of (a; 3 ) 



IV. LATTICE CALCULATIONS 

In this section, we describe salient aspects of our lattice calculation of proton matrix 
elements of the operators corresponding to moments of parton distributions. Details may 
be found in Ref. |33fl . 



A. Connected Diagrams 

Proton matrix elements of the operators in Eq. ^ are calculated by evaluating the con- 
nected and disconnected diagrams shown in Fig. |l|. Note that the term " disconnected" refers 
to diagrams in which the quarks are disconnected but of course as shown in Fig. [I], the over- 
all Feynman diagram is still connected by gluons. It is important to recognize at the outset 
that although the connected diagrams involve three propagators naively corresponding to 
the "valence" quarks in a simple hadron model, the contributions of the connected diagrams 
do not necessarily correspond to the valence parton distribution q va[ = q — q defined in 
phenomenological analyses of high energy scattering data. Rather, as shown in Eq. [I], the 
odd moments of (x n ) q and (x n )s q and even moments of (x n )& q measured in deep inelastic 
lepton scattering correspond to the sum q + q and only the even moments of (x n ) q and (x n )s q 
and odd moments of (x n )& q correspond to the difference q — q. 

It is technically much more difficult to evaluate the disconnected diagrams than connected 
diagrams and results using the eigenmode expansion technique that has recently been de- 



veloped for this purpose |34j are not yet available. Therefore, in the present work, we will 
compare the contributions of connected diagrams with flavor non-singlet combinations of ex- 
perimental results. Because the coupling of the disconnected loop to the rest of the diagram 
is flavor independent, the disconnected diagrams do not contribute to the difference between 
the moments for degenerate up and down quarks. For example, the axial charge can be calcu- 
lated directly from differences of connected contributions, qa = (Au+ Au — Ad— Ad) connecte d. 
Hence, we will subsequently compare our results with moments of flavor non-singlet combi- 
nations of the sum or differences of phenomenological quark and antiquark distributions. 

One technical difference between our calculations and those of other works is the fact 
that we use Dirichlet boundary conditions for the valence fermions in the time direction. 
With periodic or antiperiodic boundary conditions, contributions from the images of the 
sources and sinks restrict the useful source-sink separation to less than half the total time 
extent. Dirichlet boundary conditions prevent propagation from these image sources and 
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allow utilization of a larger fraction of the lattice volume for ground state operator measure- 
ments. We note that the closely related Schrodinger functional boundary conditions used in 
a calculation of the average quark momentum in the pion provide similar benefits [35] . 



B. Sources 



To facilitate measurement of proton matrix elements, we have optimized the overlap 
between a computationally efficient nucleon source and the proton ground state. The starting 
point is the interpolating field for the proton 



J 



U a C 75 d l 



u 



(12) 



where C denotes the charge conjugation matrix. This field corresponds to the non-relativistic 



it? in the non-relativistic limit. Since trun- 



quark model wave function e abc u®d\ — u^d 
cation of the lower components does not significantly reduce its overlap with the physical 
proton [35| we save an overall factor of two in computation and storage by only calculating 
propagators from the upper components of the source. 

As shown below in Section [V], if J in Eq. [T2]is a point source, the overlap with the nucleon 
ground state is of the order 10 -4 , so generalization to finite spatial extent is desirable. Since 
the source must be factorizable in order for matrix elements to be calculated from single 
quark propagators, we use distributed quark fields with adjustable spatial extent. In ref p? 



it was shown that gauge fixed sources with gaussian wave functions and gauge invariant 
smeared Wuppertal sources PB| , |3T| produced comparable results, and for computational 
convenience we use the latter. 



The smeared wave function is defined [|40] 



(N,a) 



l + aH) N ^ where H 



and ip°(n) 



J nno 



(13) 



with the spatial extent being controlled by the coefficient of the nearest neighbor hopping 
term, a, and the number of smearing steps, N. The distribution is approximately Gaussian 
with an equilibrated gauge field U producing a narrower distribution than for the free case. 
A convenient measure of the smearing is the rms radius 



/ d 3 xr 2 il>*if)~\ 2 
J d 3 xip*ip 

and figure |2] shows how r rms depends on the parameters N and a. As one expects from the 
free case and from the fact that smearing is a random walk governed by the gauge fields, the 
rms radius is approximately proportional to y/~N. Note that the size of the source is nearly 
independent of a for a > 3, at which point the constant term in Eq. [l^ becomes negligible 
relative to the hopping term. In the calculations described in section |V|, we set a = 3 and 
use N to adjust r rms to optimize the source. 

The source is optimized by maximizing the overlap between the normalized state created 
by the action of the source on the QCD vacuum = J\fJ\Q) and the normalized ground 
state of the proton |0). Denoting the momentum projected normalized eigenstates of the 
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proton by \n) and their energies by E n , the momentum projected two point correlation 
function may be expanded: 

(J(t)J(0)}=Cj:\^j\n)\ 2 e- E " t (14) 

n 

where C is an unknown normalization constant. Since one can directly measure the correla- 
tion function at zero time separation 

A = (J(0)J(0)}=Cj:\(J\n}\ 2 (15) 

n 

and reliably fit the large t behavior of the correlation function to extract the ground state 
contribution 

B = C\(J\0}\ 2 } (16) 
the probability that the source contains the proton ground state is given by 

P(0) = j = \(*j\0}\ 2 . (17) 

The importance of maximizing the overlap of the source with the proton is particularly 
clear when one considers the three point function used to calculate matrix elements: 

(J{t 3 )0{t 2 )J(ti)) =CY,(^j\n)(n\0\m)(m\ijj)e- En{t3 - t2) ~ Em ^- tl) . (18) 

n,m 

Although the contribution of a contaminant (ipj\n) ~ e, to the two point function, 
Eq. [14] , is of order e 2 , the contribution to the three point function, Eq. [18] is of order 
e {Q|0|oj - Depending on the off-diagonal elements of any particular operator, the error when 
the operator is close to one of the sources may be substantial, and the sign of the contribution 
is undetermined. Specific examples will be seen in plateau plots presented in Section [V} 



C. Sequential Propagators 

There are two alternative strategies for calculating the connected diagrams in Fig. p] with 
sequential propagators. Using a propagator from a fixed source to all points xq times the 
operator O(xo) as a new source to propagate to an arbitrary sink location allows one to 



calculate matrix elements of the operator between all source sink separations [|I3,rj,|T-H 
The matrix element is then obtained from the linear term in the time separation between 
the source and sink. The SESAM collaboration has used this variable time extent method to 



calculate the axial charge |T3[] for the configurations used in this work, and we will compare 
our results with theirs in Section [V]. However, the fact that a new set of propagators must be 
calculated for each operator makes this alternative too costly for the large set of operators 
of interest in this present work. 

The alternative we will use is to combine two forward propagators from a fixed source, 
smear them and momentum project on a specific time slice to create a sink, and use this sink 
as a new source for a backward going propagator [JPJ . The matrix element for any operator 
can then be obtained by combining forward and backward propagators with that operator. 
Because the time separation between source and sink is fixed in this approach, the effect of 
this separation has been studied and optimized as described in Section |V|. The details of 
calculating matrix elements with these sequential propagators are described in Appendix A. 
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D. Lattice Data Analysis 



To ensure that the systematic and statistical errors are well understood and controlled, 
we have used several methods to analyze the lattice measurements of matrix elements. 

In our primary analysis, which is used for all the results presented in tables and graphs, 
we define a central window for the location of the insertion of the measured operator [t source + 
St, t S i n k — St] in which contamination from excited states in the source or sink has negligible 
effect, average the ratio of three-point to two-point functions over this window, and use the 



jackknife method [46[ to estimate the average value and variance. For smeared sources, we 
used St = 3 and for tests with unsmeared sources we have also used St = 5. 

As a secondary analysis, to ensure that we have avoided systematic errors from excited 
state contaminants in the plateau region, we have also fitted the phenomenological form 

R(t) = R + J2 ( b ^ Clt + he- c ^) (19) 

i 

over the whole range between the source and the sink. The exponentials correcting for the 
contaminants are symmetric because we have used the same smearing in the source and 
sink, and the value for the matrix element is R$. 

From our data, it was only possible to determine one pair of exponentials, and we used a 
jackknife estimate on the fitted values of Rq to determine the error. For operators for which 
the errors were small enough that the exponential contaminants were well determined and 
stable, the plateau and exponential fits produced statistically consistent results with compa- 
rable errors. This ensures that our primary results do not have statistically significant bias 
from excited states. In some cases, the exponential contaminants were poorly determined 
and produced spurious fluctuations in the analysis which overestimated the actual statisti- 
cal fluctuations in the plateau region. This is the reason we selected the plateau fits as our 
primary analysis and view the exponential fits as a secondary consistency check. 

To further ensure that we understood our statistical errors, we also examined bootstrap 



distributions [47| for all our observables at the lightest SESAM quark mass. An ensemble 



of size M was created by randomly picking configurations from the original sample with 
replication, and typical distributions of values of measurements of (x)^ are shown in Fig. |3|, 
where M = 2000 and the sample sizes in the left and right figures are 25 and 204 configu- 
rations respectively. Three error bars are shown in each case: the basic jackknife error bar 
which is calculated independently of the bootstrap procedure, the 68% confidence interval 
based on the bootstrap distribution, and the statistical variance of the bootstrap measure- 
ments. The three arrows correspond to the jackknife result, the median of the distribution, 
and the mean of the distribution. This figure shows that for a sample size as small as 25, 
the distribution is far from Gaussian, the three error bars differ significantly, and even the 
arrows disagree slightly. For the large sample size 204 used in our full data analysis, every- 
thing is consistent: the distribution is nearly Gaussian, the three errors are consistent, and 
the three arrows are equivalent. Thus, the bootstrap analysis gives added assurance that 
the statistics are understood and under control. 
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V. LATTICE RESULTS 



This section presents the results of our calculations of moments of quark distribution 
functions in full and quenched QCD, as well as relevant tests of the lattice technology and 
consistency checks. 



A. Lattice Gluon Configurations 

The parameters specifying the full QCD, quenched, and cooled configurations used in 
our calculations are tabulated in Table |. The full QCD configurations were calculated 



using the hybrid Monte Carlo algorithm by the SESAM collaboration []16[ at j3 = 5.6 and 
by the SCRI collaboration [^] at (3 = 5.5. The quenched configurations were calculated 
at MIT at (3 — 6.0 to be directly comparable in lattice spacing to the unquenched (3 = 5.6 
configurations, and the cooled configurations were obtained by applying 50 cooling sweeps 
to the SESAM configurations as described later. 

The practical limitations associated with these configurations pose significant problems 
in extrapolating to both the chiral and continuum limits as required to compare with phe- 
nomenological data. 
Chiral extrapolation 

A common criterion to keep the pion correlation length significantly smaller than the 
physical size of the lattice and to avoid unphysical interactions with periodic images is to 
require that the pion Compton wavelength be less than one-fourth the spatial dimension 
of the lattice. Because resources limit us to 16 3 x 32 lattices, at /3 — 5.6 where the lattice 
spacing a = 0.091 fm, the lattice dimension is only 1.46 fm so that the pion mass must be 
greater than 540 MeV. However even for the lowest SESAM quark mass, which exceeds this 
criterion, there are significant discrepancies between mass measurements on 16 3 x 32 and 
24 3 x 48 lattices |4~9|, indicating serious finite volume effects. Hence, in our final analysis, 



we have only included the heaviest three SESAM quark masses in our results. 

As discussed in detail in the final section, much of the physics of the pion cloud of the 
nucleon is omitted on a lattice of size 1.46 fm with pions heavier than 540 MeV, so it is 
unreasonable to expect naive linear extrapolation of calculations in this regime to accurately 
include the quantitatively important physics of the pion cloud. However, since we do not 
presently have data in a regime in which we can fully determine the parameters of an 
extrapolation incorporating the behavior known from chiral perturbation theory, most of 
our data analysis will be based on naive linear extrapolation. This extrapolation is still 
useful because it crisply frames the issues that need to be confronted in comparison of ab 
initio lattice calculations with experiment. 
Continuum extrapolation 

Due to limitations in presently available configurations for Wilson fermions, extrapolation 
to the continuum limit is even more problematic. Table |y] summarizes parameters and 
lattice spacings for configurations that are presently available for dynamical Wilson fermions. 
Since we are calculating nucleon properties, we believe it is most consistent to set the scale 
using the lattice spacing determined from the nucleon mass, a at- Unfortunately, there is a 
large discrepancy in determinations of a at, and for qualitative purposes we have assumed 
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the behavior given by the LANL results, since this is the only data set in which one group 
calculated consistently at three different (3. 

One goal has been to use a comparison of calculations of moments of quark distributions 
using SCRI configurations at (3 — 5.5 and SESAM configurations at (3 = 5.6 to obtain an 
indication of the finite lattice size dependence, and we report results from both /3's below. 
However, the problem of finite size effects is shown in Fig. f|, where we present scaled values 
of m 2 at different values of (3 for the data sets in Table [TV]. The scaled mass squared, 
(am^) 2 x (a 5 . 6 /a) 2 , converts the mass in different lattice units at each a to fixed lattice units 
at [3 — 5.6 for ease of comparison. Thus, the ordinate is proportional to m 2 in physical 
units, and 0.05 corresponds to m 2 = 0.24 GeV 2 or m n = 490 MeV. The dotted line shows 
the point at which the pion Compton wavelength equals one fourth of the spatial dimension. 
Given that the lowest SESAM point has significant finite size errors, we expect the lowest 
two SCRI points to have comparable contamination from finite size effects. Since there are 
no SCRI points comparable to the highest SESAM points, there is no way to delete the 
lowest SCRI points and perform an extrapolation comparable to the SESAM extrapolation. 
Hence, we believe that the differences between moments of quark distributions calculated 
with the SCRI and SESAM configurations is an undetermined combination of lattice spacing 
errors and finite volume effects. New dynamical Wilson quark calculations are underway at 
(3 = 5.5 and (3 = 5.3 to separate these finite volume and finite lattice spacing effects. 



B. Source optimization 

To obtain as much physics as possible from the SESAM and SCRI configurations, the 
nucleon sources have been optimized in their overlap with the ground state and with respect 
to their separation in Euclidean time. 

As discussed in section [IV B| , the rms radius, (r 2 )?, characterizing the spatial extent of 



a Wuppertal source was varied to maximize the overlap of the source with the ground state 
proton wave function, P(0) = |(\l/j|0)| 2 . Figure [5] shows the overlap, with jackknife errors, 
between the smeared source and the proton ground state as a function of the source rms 
radius in lattice units, calculated at (3 = 5.6, n = 0.1575, and smearing parameter a = 3. 
The number of smearing steps, N, to produce the rms radii in this plot ranges from 20 at 
the lowest non-zero point to 250 for the highest. It is interesting that the overlap increases 
from 6 x 10~ 5 for a point source to 0.5 for an rms radius of the order of 4.5 lattice units. 
Clearly, this four-order-of-magnitude increase is a dramatic aid in improving the plateau 
for measurements. Based on these and analogous calculations at — 5.5, the optimal 
Wuppertal smearing occurs at a physical rms radius of 0.4 fm, which corresponds to a= 3 
and N = 50 at (3 = 5.6 and a= 3 and N = 30 at (3 = 5.5. 

Note that this overlap calculation opens the possibility of performing instructive varia- 
tional studies of the nucleon ground state. Already, a very simple quark and gluon wave 
function made up a product of three gaussian single-particle wave functions smeared with 
surrounding links of glue has a 50% overlap with the full ground state and is measured 
with good statistical accuracy. Hence, it should be straightforward and practical to inves- 
tigate variationally the extent to which the wave function can be improved by the addition 
of changes in the quark wave function, such as the addition of diquark correlations, or of 
changes in the gluon wave function. 
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The optimal separation between source and sink is a compromise between two compet- 
ing effects. As the separation increases, the size of the usable plateau increases, reducing 
systematic error and in principle allowing for increased statistical accuracy from averaging 
over an increasingly large region. However, the propagators themselves are exponentially de- 
creasing with separation and the fractional error with which they can be measured increases 
exponentially. (This is analogous to the exponential error growth in the measurement of 
Wilson loops or Polyakov lines, but because the measurement of the three-point function is 
intrinsically non-local, it is not amenable to the error reduction techniques of Liischer and 
Weisz g). 

To determine the optimal compromise between these two competing effects more quan- 
titatively, plateaus were calculated as a function of the source-sink separation for point 
sources, and the results for (x)^ using 100 configurations at /3 = 5.6 and k = .1575 are 
shown in Fig. |] for time separations AT = 12 and 14. The errors for AT = 14 are 50% 
larger than for AT = 12, and become prohibitive for larger separations. This is consistent 
with the results of ref |H|, that for a 16 3 x 32 lattice in quenched QCD with (3 = 6.0, which 
has a comparable lattice spacing, the error bars grow significantly for AT > 13. Although 
the plateau region for AT = 14 is hardly present at AT = 12 for point sources, one may 
still extract the correct matrix element using the exponential fit method described in section 
LV D|. Furthermore, using optimal smeared sources to enhance the ground state overlap by 



four orders of magnitude decreases the excited state contaminants to the point that a well 
defined central plateau is recognized for AT = 12 as shown below. Hence, to reduce the 
overall statistical errors in our calculations, it is optimal to use a physical separation of 
approximately 1.1 fm which corresponds to AT = 12 for SESAM configurations at (3 = 5.6 
and to AT = 10 for SCRI configurations at j3 = 5.5. 

The final quality of the plateaus used for measurements of operators in a zero mo- 
mentum ground state are shown in Fig. ^|, where we plot measurements of the operators 
(1)a 9 , (^)aL an d (1)<5 9 as a function of Euclidean time for an ensemble of 200 SESAM 
configurations with (3 = 5.6 and n = 0.1560. The improvement in the plateau for (x)^ 
produced by smearing is seen by comparing the right panel of Fig. ^] and the upper left 
panel of Fig. [5]. Note also, as emphasized in connection with Eq. that the sign of the 
exponential contaminant near the source or sink may be either positive or negative, with 
a small negative contribution arising in the case of the axial charge, (1)a<j shown in the 
upper right panel of Fig. 0. These well-defined plateaus are typical of all our measurements 
of operators in a zero momentum ground state and show why we obtain consistent results 
with plateau and exponential fits and satisfactory statistics. In the case of projection onto 
non-zero momentum ground states, errors are substantially larger. Although we also report 
these non-zero momentum results in this work, the errors are generally unsatisfactory and 
additional calculations will be carried out to improve their statistics. 

One other optimization was investigated, but turned out to produce minimal im- 
provement. There is considerable freedom in selecting lattice operators from the irre- 
ducible representations tabulated in Table 0. For example, instead of the basis vector 
qj^D^q — 3(971 Diq + 972-D29 + q^izDzq) of representation 3^ , one could equally well chose 

q^iDiq — q"fiDiq for i — 1, 2, or 3. Since the direction i = 3 is inequivalent to the other two 
spatial directions because the spin is polarized in the 3 direction, one can find numerically 
the optimal combination of the basis vectors that minimizes the variance in the measure- 
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ment of (x)^. However, in the end, the reduction in the statistical error is only about 7%, 
so this refinement was not implemented in our production calculations. 

C. Consistency checks 

This section summarizes a number of consistency checks that have been performed to 
ensure the reliability and accuracy of our calculations of moments of quark distributions. 
Sources 

One potential pitfall of optimizing sources to increase the size of the plateau region is 
the possibility that although excited state contaminants appear to have been diminished, 
unobserved systematic errors may have sneaked into the plateau measurements. To ensure 
that such systematic errors have not been introduced into our present calculations, we have 
evaluated the operators of interest for various combinations of source and sink combinations. 
Figure § compares the results for (x)f\ (1)a<3, and (1)^ calculated using point-point, point- 
smeared, smeared-point and smeared-smeared source sink combinations, showing that in 
all cases the results are statistically consistent. Forty configurations at (3 = 5.6 were used 
with a source-sink separation of 12, and the smeared sources had N = 20 smearing steps. 
Each cluster of four error bars corresponds to four different window sizes over which the 
measurements were averaged, characterized by the number of points omitted at the source 
and sink. The case shown in Fig. [F], denoted (3,3), corresponds to omission of 3 points 
at the source and 3 points at the sink and is shown at the left of each cluster. With a 
separation of 12 lattice spacings, there are 13 points in all, of which 6 are deleted leaving 
a measurement plateau of 7. The remaining error bars in each cluster correspond to (3,5), 
(5,3), and (5,5) points omitted at the source and sink respectively. Thin lines denote cases 
in which the window comes too close to either a point sink or point source and overlaps with 
a region with contamination from excited states and the error bars denote jackknife errors. 

Figure |9| shows a related comparison of the measurements of the same three observables 
in which the source and sink are smeared with N = 0, N = 20, and N = 100 smearing 
steps. The window sizes within each cluster of error bars are the same as in the previous 
case. Again, the observables are consistent for all three cases, although the error bars become 
substantially larger for the largest smearing because the size of the paths of link variables 
generated by the smearing increases significantly. 
Boundary conditions 

The Dirichlet boundary conditions in Euclidean time we have used enable us to calculate 
the exponential decay of two point functions far beyond the midpoint of the lattice without 
the usual contributions of the propagation of the parity partner of the nucleon from first 
images in the time direction. Only within a few lattice points of the edge of the lattice 
do artifacts associated with reflection from the boundary (or equivalently, negative image 
charges) become significant. To verify that the choice of boundary condition does not affect 
the measured values of hadronic matrix elements, Fig. [Tl] shows the results of measuring the 
operators (x)^\ (I)a,?, (1)<5<?; an d (x)A q using Dirichlet and periodic boundary conditions. 
Note that although the results are statistically consistent, our choice of Dirichlet bound- 
ary conditions significantly reduces the error bars for some operators, such as (x)^' and 

(1}a<7, (Im- 
precision 
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In our production calculations, the conjugate gradient inversion to calculate quark propa- 
gators was performed in single precision (to increase performance on the ES-40 cluster) using 
a perturbative stopping residue r^ in = 10~ 14 and four restarts |HJ. To verify that single 
precision with the restart procedure provided adequate precision, we show the discrepancy 
between calculations of the axial charge, (1)aij, calculated in single and double precision as a 



function of the residual r^ in in Fig. O. We note that all discrepancies are already negligible 



at residue r 2 min = 10~ 12 , providing a two order of magnitude safety margin. 
Autocorrelation functions 

Autocorrelation functions were calculated for lattice measurements of (x)}® and (x)& q 
using the SESAM ensemble at k — 0.1575 consisting of 200 configurations separated by 
25 Hybrid Monte-Carlo trajectories. Correlations for all non-zero configuration separations 
were statistically consistent with zero, and all indications suggest that the 200 configurations 
are indeed independent for all the degrees of freedom relevant for the present measurements. 

In contrast to the SESAM configurations, the SCRI configurations are only separated by 
10 trajectories. To minimize the possible effect of correlations, the physical location of the 
source and sink was shifted by 16 time slices on all odd configurations relative to those on 
even configurations. Autocorrelation functions for 200 SCRI configurations at k — 0.1600 
yielded a result statistically inconsistent with zero only for separation by two configurations, 
in which case the correlation was approximately 0.2 ± 0.15. Thus for the degrees of freedom 
relevant to our proton matrix elements, configurations separated by 10 Hybrid Monte-Carlo 
trajectories with sources and sinks displaced by 16 time slices are independent, whereas those 
separated by 20 Hybrid Monte- Carlo trajectories with the sources and sinks at the same 
location still have non-vanishing correlations. The relative independence of measurements 
at different source and sink positions provides an attractive means of extracting additional 
physics from the available full QCD configurations that will be exploited in subsequent 
calculations. 

Two and three point functions 

The nucleon two point function may be calculated in two ways: by combining at the 
sink three forward propagators calculated from the source, or by combining a backward 
propagator with the sequential source. The agreement within roundoff error of these two 
results was used not only as a check on source methodology, but also at run time as a test 
to verify that stored forward and backward propagators had not been corrupted. 

Another check is to verify that the spatial integral of the matrix element of the con- 
served nonlocal vector current (the Noether current for the discrete Wilson lattice action), 
(P\ J d 3 x Jq' NL \P) I (P\P), equals 2 or 1 for up and down quarks respectively. This relation 
must be satisfied configuration by configuration and is also used as a run-time test of the 
configurations. 

An additional test of the source and sink construction is to project the sink onto a fixed 
position in the sink time-slice, rather than momentum projecting it, and to replace all link 
variables by their mirror images. We have verified that the resulting three-point function is 
the mirror image of the original result configuration by configuration. 

Finally, as discussed previously, matrix elements of operators may either be calculated 
with the variable time extent method, in which one sums over all times at which the operator 
insertion could occur and determines the slope of the dependence on t S i n k, or by measuring 
the plateau using the fixed source-sink separation as done in this work. A consistency check 
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that is satisfied configuration by configuration is that the sum of the three-point function over 
the intermediate time at which the operator is measured for a fixed source-sink separation 
should be equal to the variable time extent three point function for the fixed source-sink 
separation. In addition, we have verified that our measurements of the axial charge using 
the fixed source-sink separation method for the three largest values of k agree quantitatively 
with those calculated by the Wuppertal group using the variable time extent method |14[ . 



On the basis of these and other consistency checks, we are confident that the calculations 
of moments of quark distributions reported in this work are both reliable and accurate. 



D. Quenched results 

In order to investigate the effects of quenching quantitatively, it is essential to compare 
quenched and unquenched calculations at equal lattice spacing using precisely the same 
computational methodology. Hence, we have calculated moments of polarized and unpo- 
larized parton distributions in quenched lattice QCD at (5 — 6.0, which is comparable to 
full QCD at (3 = 5.6. These calculations also afford the opportunity to check our results 



with the extensive calculations by the QCDSF collaboration |[lC||lq , |52| - ^1 . For Wilson 
fermions at (3 = 6.0, the critical k is k c = 0.1572, and we have performed our calculations at 
k = 0.1550, 0.1540, 0.1530. Two of these values of k coincide with values used by QCDSF, 
facilitating direct comparison of unextrapolated as well as extrapolated results. In addition, 
the pion masses corresponding to these values of k overlap the range of pion masses for the 
SESAM configurations, providing a meaningful comparison of quenched and unquenched 
QCD. 

Linear extrapolations of the quenched moments of unpolarized up and down quark dis- 
tributions (x))p and (x 2 ) q are shown in Figs. [12] and [13]. Corresponding results for the 



quenched moments of polarized quark distributions (1)aq, and (x)^ are shown in Figs. |14 
and [15]. The complete numerical data for these and other moments of quark distributions 
are tabulated in Table |V|. All results are presented in the MS scheme with fx 2 = 4 GeV 2 
and extrapolated linearly with a least squares fit to the chiral limit. In these figures, the 
measured values of m 2 at each k were converted to physical units using - = 2 GeV. 

Our results for all operators agree within statistics with those of QCDSF, and Figures 



12] ^1~5| show typical comparisons. Note that for the three operators that can be evaluated in 
a zero-momentum proton ground state, (x)f\ (1)a<3, and (x)^ q the statistical errors in the 
extrapolated results are small and allow for meaningful comparison with full QCD results 
and with phenomenology. In contrast, the errors associated with the operator (x 2 ) q , which 
requires projection onto a state with non-zero momentum, are so large that meaningful com- 
parison will require higher statistics. Since QCDSF used on the order of 1000 configurations 
whereas we only used 200 configurations, it is consistent that our error bars on individual 
points are twice as large as theirs for this operator. The fact that the QCDSF error bars 
are not half the size of ours for three matrix elements calculated in a zero momentum state 
as well is explained by the fact that the errors they reported for these matrix elements 
are somewhat larger than jackknife errors because their original analysis used a full error 
correlation matrix from which low eigenmodes were removed ||56||. 
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E. Full QCD Results 



A central result of this work is the fact that for the range of quark masses accessible in 
the present calculations, quenched and full QCD calculations at comparable lattice spacing 
agree within statistical errors. The moments of quark distributions in full QCD at (3 — 5.6 



calculated using SESAM configurations are shown in Table VI. Note that to avoid finite 



volume effects, moments at the lightest quark mass are not included in the extrapolation. 
Detailed comparison with Table [V] shows that in the case of all the operators we have cal- 
culated, the full QCD results at (3 = 5.6 are statistically consistent with those calculated 
in quenched QCD at (3 = 6.0. Furthermore, in the case of operators that can be evalu- 
ated in zero-momentum nucleon states, the statistical errors are sufficiently small that any 
differences between the linearly extrapolated quenched and unquenched results are small 
compared with the discrepancies with experiment discussed below. The level of agreement 
between quenched and full QCD of (x 2 ) q , (1)a^ 5 ( X /A q i an d (l)<s 9 is also shown in Figs [TB|, 



T7L [18], and [19] respectively. Hence, the present work strongly rules out the conjecture that 



the serious discrepancies between quenched calculations and phenomenology arose from the 
omission of the effect of dynamical quarks. As noted previously, the magnitude of the errors 
in observables measured in nucleons with non-zero momentum is qualitatively larger, and no 
strong conclusions can be reached for these observables until higher statistics calculations 
are carried out. 

In addition to the calculations at (3 — 5.6, which correspond to the same lattice spacing 
as our quenched calculations, we have also calculated the same operators at (3 = 5.5 using 
the SCRI configurations in an attempt to study the approach to the continuum limit. The 
results are tabulated in Table |V11| . Note that by Eq. 10, the shifts in a and j3 in going from 
(3 = 5.6 to (3 = 5.5 produce changes in the renormalization constants that are negligible on 
the scale of the present statistical errors, so these changes in renormalization factors have 
been omitted. Unfortunately, as pointed out in connection with Fig. |], the two lightest quark 
masses lie in a regime in which we expect significant errors due to finite- volume effects, so the 
chiral extrapolations are physically suspect. Even so, the bulk of the operators evaluated 
at (3 — 5.5 and (3 = 5.6 agree within statistics. The only exceptions, which are only 
slightly beyond one standard deviation, are (x 2 ) u , (1)a«; ( x )au > (1)<5«> an d ( x )sd- Because 
of the uncertainty in finite volume effects and small lever arm, we have not attempted to 
extrapolate the present data to the continuum limit. Higher statistics calculations at smaller 
quark masses for several coupling constants are presently underway to more seriously address 
the behavior in the continuum limit. 



F. Cooled Results 

A fruitful strategy to obtain insight into hadronic physics is to use the numerical evalu- 
ation of the sum over all quark and gluon configurations contributing to the path integral 
to isolate those paths that dominate the action. In recent years, this approach has provided 
strong evidence that in QCD with light quarks, topological excitations of the gluon field, 
which in the semiclassical limit correspond to instantons, play a major role in hadron struc- 
ture. By minimizing the action locally in a process known as cooling []5?|], the instanton 



content of the quenched and full QCD vacuum has been extracted [^8]-|60||. Comparison of 
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hadronic observables calculated with all gluons and those obtained using only the instantons 
remaining after cooling has demonstrated qualitative agreement for hadron masses, quark 
distributions, and vacuum correlation functions of hadron currents |JT] . Calculation of the 
lowest quark eigenmodes has revealed zero modes correlated spatially with the instantons 
and truncation of the quark propagators to the zero mode zone has produced the full strength 
of the p and 7r contributions to vacuum correlation functions |62|j63| . 

Hence, to obtain further insight into the structure of the proton, we have used cooling 
to remove essentially all the gluon fluctuations except for instantons from the full QCD 
configurations, and compared the quark distributions calculated in full QCD and including 
only the contributions of instantons. The role of instantons is particularly interesting in 
considering the spin structure of the proton, since the 't Hooft instanton interaction is 
the only vertex in QCD that directly removes helicity from valence quarks and transfers 
it to gluons and quark-antiquark pairs and is therefore a natural mechanism to explain 
the so-called "spin crisis" . Our calculation of the contributions of instantons to 

matrix elements of operators related by the operator product expansion to deep inelastic 
scattering also closely parallels the direct calculation of instanton contributions to deep 
inelastic scattering in refs. [|67| |TQ|1 . The present work extends and substantiates earlier 
exploratory investigations [I2|. 

To remove the essential non-instanton related gluon fluctuations without producing un- 
necessary annihilation of instanton-antiinstanton pairs, we have cooled the full QCD SESAM 
configurations at (5 — 5.6 using 50 cooling steps. This amount of cooling corresponds to 
roughly 25 cooling steps for (3 = 5.7, which was an effective amount of cooling in ref flBT 



Furthermore, calculations of nucleon two-point functions using the SESAM configurations 
have shown that the results with 25 and 50 cooling steps differ negligibly [[51 . 

Because of the smoothness of cooled configurations, small statistical errors are obtained 
using an ensemble of 100 configurations, and previous investigations have shown that the 
chiral dependence on quark mass is quite linear in the region of interest. Hence, we cooled 
100 configurations at K sea = 0.1570 and at K sea = 0.1560 for our comparison with full QCD. 
To make the cooled chiral extrapolation comparable to the uncooled case, for each n sea , 
we selected the cooled valence quark mass, n va i, such that the ratio was the same as 
the ratio in full QCD. The resulting cooled values are K va j, = 0.1246 for n sea = 0.1570 , 
where = 0.523, and K vcd = 0.1235 for n sea = 0.1560, where = 0.476, and the full 
QCD masses were taken from ref |16[]. Note that, because we are comparing two theoretical 
calculations in the same physical volume, we have not discarded the lightest quark mass case 
for full QCD out of concern for finite volume effects, but rather have included it to enable 
comparison in the regime of the lightest quark masses where dominance by the zero modes 
associated with instantons should be most pronounced. 

Since the high frequency quantum fluctuations are removed by cooling, we set all the 
renormalization constants, Z, to one in order to compare cooled results with the full QCD 
results. Although we are aware of no rigorous argument as to formulate renormalization in 
the presence of cooling, this approximation appears the most physical and, for example, is 
quite sensible for tadpoles, where we expect Uq = (^TrU) 1 ^ ~ 1. 

The moments of quark distributions calculated in cooled configurations are tabulated in 



Table |VIII| . In addition, extrapolations in full and cooled QCD for the operators that can 
be calculated with high statistics in a nucleon state at zero momentum, (x)^\ (1)a 9 , (x)aI, 
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(l)sq, and d\, are compared in Figs. ^T|, [22|, |23|, and [24| respectively. 

It is striking that the extrapolated cooled and full QCD results agree so closely for all 
the twist two operators that correspond to moments of quark distributions. This detailed 
agreement, generally within error bars but always within two standard deviations, provides 
strong support for the physical picture that the propagation of light quarks in the nucleon 
is strongly dominated by the instanton content of the gluon configurations. It is interesting 
that in the cases of (x)^ and (x)^ q , the cooled and uncooled results agree best in the region 
of light quark masses, where we expect zero mode dominance to be most pronounced, and 
differ much more significantly at heavy quark masses where they have no reason to agree in 
detail. 

Having observed consistency between cooled and full QCD results in the cases above 
where one expects agreement on the basis of instanton physics, it is also interesting that 
cooled and full QCD results differ by an order of magnitude for the twist-three operators 
d\ and d 2 where we expect them to disagree dramatically because of operator mixing. Re- 
call that in connection with Eq we pointed out that for Wilson fermions, the operator 

Ts7[o-I?{Ati] • ■ ■ Dti n } mixes with the lower dimension operator ^757[ct7{ w ] • ■ ■ D^}- Hence, one 
expects that our lattice measurements of d\ and d 2 are contaminated to a large extent by 
operator mixing, and indeed a calculation of the non-perturbative mixing |J3J yields a large 
change in the extracted value of d 2 . As argued above, cooling removes the short wavelength 
fluctuations responsible for renormalization and mixing, so as a result we would expect the 
mixing to be reduced essentially to zero. This is precisely what is observed in Fig. 53 and in 
Table |VIII| where the cooled measurements extrapolate to values close to zero and the full 
QCD measurements are an order of magnitude larger. 



G. Comparison with Phenomenology 

This section presents the evidence that linear chiral extrapolation of full lattice QCD 
results in the regime of quark masses accessible in this work is seriously inconsistent with 
experimental measurements of light cone quark distributions in the nucleon. To make the 
argument as strong as possible, it is useful to first compare our calculations with three 
other lattice calculations. Hence, the quenched and full QCD calculations of this work are 
compared with other related lattice calculations and with phenomenology in Table ITX. 



We have already commented on the comparison of our quenched calculations with those 
of the QCDSF collaboration, and their results from refs pH| , |52| -|55B are tabulated in the first 
column to be compared with our quenched results in the fourth column. Here, complement- 



ing the plots in Figs [12] - [15], one sees detailed agreement of results calculated at the same 
j3 and k, strongly supporting the accuracy and consistency of both calculations. 

The second column shows the extrapolation to the continuum limit by the QCDSF 
collaboration of several moments |TJ[]. The agreement of these extrapolations with the first 
column clearly shows that in the case of quenched QCD, finite lattice spacing effects are 
small compared with the discrepancies with experiment. Given the close agreement between 
full QCD and quenched calculations in this regime of quark masses, there is no reason to 
believe that finite lattice spacing effects are substantially larger in our full QCD calculations 
than in the quenched case, so it would be unreasonable to attribute discrepancies with 
experiment to finite lattice spacing effects. 
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The third column shows the contributions to the axial charge by the SESAM collabora- 
tion using the same gluon configurations but completely different measurement technology 
TM to be compared with our results in the fifth column. As in the case of the quenched com- 



parisons, we believe this provides additional confirmation of the accuracy and consistency 
of the present measurements. 

The final column shows the moments that were calculated numerically from phenomeno- 
logical fits to the world supply of deep inelastic lepton scattering data and other high energy 
scattering data. Data is conveniently accessible on the web from each of the major collab- 
orations, and we have calculated unpolarized moments using each of the unpolarized data 
sets, CTEQ U, GRV [0, and MRS [Q] and calculated polarized moments using both of the 
sets GRSV 5U, and GS ||. Unfortunately, it is difficult to provide quantitative estimates of 
systematic or statistical errors. Experimental data is available only over part of the neces- 
sary range of momentum fraction, x, and model assumptions are invoked to parameterize 
parton distributions at large x consistent with known sum rules and physical constraints. 
In addition, no error correlation matrix is provided, so it is not even possible to calculate 
the error in a moment that arises from the known statistical errors in the measurements. In 
order to get some indication of the possible errors, we have calculated the moments using 
each of the unpolarized or polarized data sets. In the table, we tabulate the average value 
and give the maximum difference between values for all the relevant data sets in parentheses. 
Note that these differences are small compared to the statistical errors in the corresponding 
lattice measurements, so we believe the phenomenological uncertainties are small compared 
to the discrepancies with lattice extrapolations discussed below. Since, as argued previously, 
the disconnected diagrams do not contribute to the flavor non-singlet combination {0) u -d, 
we have tabulated the differences between the up and down quark contributions for the 
relevant combination of q ± q for comparison with the lattice calculations. 

Table |X| reveals a clear discrepancy between linearly extrapolated lattice calculations 
and phenomenology. The momentum fraction (x) q is a fundamental property specifying 
the fraction of the total momentum carried by a particular quark flavor. The non-singlet 
momentum fraction is overestimated by more that 50%, with full QCD yielding 0.25 - 
0.29 compared with the phenomenological result 0.15. The second and third moments are 
similarly overestimated by linear extrapolation. 

Another important quantity is the nucleon axial charge, (1)au-a<2 governing (3 decay. 
Here again, one finds major discrepancies with linear chiral extrapolation, with quenched 
or unquenched calculations yielding 1.0 - 1.15 compared with the experimental value 1.26, 
corresponding to a discrepancy of 10 to 25%. 

Since these and other comparisons in Table [TX] show a clear discrepancy between phe- 
nomenology and linear chiral extrapolations and because of the evidence summarized pre- 
viously for the consistency and accuracy of the lattice calculations at heavy quark masses, 
we believe that the fault lies with the linear chiral extrapolation as discussed in the next 
section. 



H. Chiral Extrapolation 

Superficially, one could imagine that since deep inelastic scattering involves large Q 2 , it 
might be dominated by short distance behavior of the nucleon wave function and not be 
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strongly influenced by the long distance pion cloud. However, as we have seen, the operator 
product expansion relates moments of structure functions to a tower of local operators to be 
evaluated in the nucleon ground state. Physically it is clear that the pion cloud should play 
a major role in the nucleon matrix elements of these operators, and that these contributions 
are strongly suppressed by the heavy quark masses and small spatial volumes in which we 
have been forced to work. Whereas there is clearly a heavy quark regime in which matrix 
elements vary nearly linearly with the quark mass or m 2 , the behavior becomes highly 
non-linear at quark masses sufficiently light that a substantial pion cloud is produced. 

Complementary to the linear heavy quark regime, there is a second regime near the chiral 
limit in which the behavior is again simple and is specified by chiral perturbation theory. 
Here, the physics is described by an effective chiral theory based on the would-be Goldstone 
bosons, and the leading non-analytic behavior can be calculated uniquely. For example, of 
direct relevance to this work, the leading non- analytic behavior of (x n ) u ^d is |7l| , [72[| 

{X n )u- d ~ On [l ~ {3 f^ f l]T l ^] ■ (2°) 

Unfortunately, at present, there is no full analytical theory joining the chiral regime 
with the heavy quark regime of the present lattice calculation. Hence to explore the chi- 
ral extrapolation of lattice QCD calculations to the physical region, it is useful to use a 
physically motivated extrapolation formula incorporating the correct behavior in the chiral 
limit. We therefore fit the lattice data with the extrapolation formula of Ref. ||18|| , in which 



a phenomenological cutoff /i is introduced in the nonanalytic term to specify the size of the 
source generating the pion cloud and the usual analytic term in m\ is included: 

<*">.-, ~ ^ - ^ ln (<f^)l +K< (21) 

Physically, it is reasonable that momenta in the pion loop should not become infinitely 
large as they would in the presence of a point source, but rather be cut off at a scale charac- 
teristic of the size of the valence quark core of the nucleon, and alternative parameterizations 
in terms of a form factor give equivalent results. 

The result of using this extrapolation formula for the difference between the up and down 
quark momentum fraction (x n ) u -d is shown in Fig. Although this extrapolation is only 
valid in full QCD, since full QCD and quenched results are equivalent in the regime of our 
calculations, to improve statistics, we have also included our quenched results. The heavy 
solid curve is the result of a least-squares fit of ai and bi in Eq. |21] with fixed \i = 550MeV 
and the the light solid lines indicate the jackknife error band. Here we see clearly that 
the extrapolation formula containing the leading chiral behavior is consistent both with 
the lattice measurements and the experimental data. The fact that the effect of the chiral 
logarithm is much larger in this matrix element than in more familiar mass measurements 
is again a manifestation of the variational theorem. Making an order e error in the pion 
cloud of the wave function makes an order e error in matrix elements of general operators 
but only an order e 2 error in the mass. The obvious problem with the present argument is 
that although a cutoff of order 0.3 fm is physically reasonable, it has not been calculated 
from first principles and may be regarded as a single free parameter introduced to fit a single 
experimental measurement. 
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To see that this single cutoff parameter resolves the discrepancy with experiment for a 
number of observables, the results of using Eq. 21 for the first three moments of the differ- 
ence between the up and down quark density distributions {x n ) u _d are shown in figure |26|. 
Because of the larger error bars in the higher moments, we have also included the quenched 
QCDSF results ]TD| to provide a larger lever arm in the least squares fit. As seen in the 
figure, the single value of \i provides simultaneous agreement with experiment for all three 
moments. As shown in Ref. []18|| , this strong chiral behavior of the three lowest moments is 
also observed in chiral bag models. Furthermore, similar extrapolation with a comparable 
cutoff also reconciles the strong discrepancy between linear extrapolation of lattice results 
and experimental nucleon magnetic moments JO]. Thus, a consistent picture is emerging 



concerning the importance of large effects of the pion cloud in chiral extrapolations. Unfor- 
tunately, because spin-dependent structure functions also involve significant contributions 
from Delta excitations and chiral perturbation theory becomes less well controlled, we do 
not presently have a corresponding physical interpolation formula for (x n )Au-Ad- 

The curves in Fig. ^ suggest that precision measurements down to m% = 0.05 GeV 2 are 
required to determine the parameters of the chiral extrapolation and thereby provide reliable 
extrapolation of moments of quark distributions. The computational resources required for 
such calculations may be estimated using the cost function J73[] obtained by the SESAM 
collaboration, which provided the gluon configurations. For present purposes, the number 
of floating point operations per independent gluon configuration, N in Teraflops-years, may 
be conveniently written 



N ~ .038 



4.55 r Qg-i 7.25 



.3 



2.7 



Because the spatial derivatives and non-zero momentum projections required to calculate 
moments of structure functions require high Monte Carlo statistics, it is necessary to calcu- 
late of the order of 400 independent configurations. Including equilibration and calculation 
at higher quark masses, the total computer time is approximately twice that required for 
400 configurations at the lowest quark mass. Hence, a calculation with a lattice spacing 
a = 0.1 fm and — = 0.3 would require approximately 8 Teraflops-years, i.e., dedicated use 
of a computer that sustains 8 Teraflops on QCD for one year. Such resources should become 
available in the next generation of lattice QCD computers. 



VI. SUMMARY AND CONCLUSIONS 

This work has presented the first calculation of the moments of light cone quark distribu- 
tions in full QCD. The methodology has been presented in detail and validated by numerous 
consistency checks and comparison with other relevant lattice calculations. 

One major result of this work is the close agreement of full QCD and quenched calcula- 
tions for quark masses corresponding to above 500 MeV. This agreement rules out the 
conjecture that discrepancies with experiment could be attributed to quenching effects. 

A second salient result is clear evidence that linear extrapolation of full QCD lattice 
results from pion masses above 500 MeV is in serious disagreement with experiment, ranging 
from the order of 50% for the quark momentum fraction to 10-25% for the axial charge. 
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Motivated by the inadequacy of linear chiral extrapolation, we have shown that extrap- 
olation including the leading nonanalytic behavior of chiral perturbation theory has the 
potential to yield results consistent with phenomenology. We have explicitly shown that 
an extrapolation formula with a single phenomenological cutoff simultaneously fits the first 
three moments of the quark momentum fraction. 

Finally, we have shown qualitative agreement between full QCD and cooled lattice con- 
figurations retaining only the contributions of instantons, providing additional evidence for 
the role of instantons in light hadron structure and of zero-mode dominance. 

This work points the way for a number of promising steps in our continued quest to 
understand hadron structure from first principles using lattice QCD. One should clearly 
undertake a systematic program of using a partially quenched chiral expansion for extrap- 
olation and measure the relevant parameters of the effective chiral theory. To complement 
the flavor non-singlet matrix elements of this work, we need to calculate the disconnected 
diagrams required to compare flavor singlet matrix elements with phenomenology. The con- 
tinuum limit of the SESAM results should be explored by supplementing the existing SCRI 
configurations with additional quark masses so that calculations at f3 — 5.6, 5.5, and 5.3 may 
be compared and extrapolated. It is desirable to undertake nonperturbative renormalization 
without Gribov ambiguities associated with gauge fixing, and to this end we note that the 
Schrodinger functional method has now been used to calculate the average quark momen- 



tum in the pion p5] , |76T [78|1 . Quenched calculations with chiral fermions in a spatial volume 
of 3.2 fm should be carried out to extend the pion mass down to 250 MeV and to remove 
the problem of operator mixing in the twist-3 matrix elements. Finite volume formulae for 
operator matrix elements should be derived to correct residual finite volume effects. As 
lattice results come closer to phenomenology, it will be important to have quantitative un- 
derstanding of the errors in moments of structure functions arising from systematic errors, 
as well as the error correlation matrix describing statistical errors. Finally, with the advent 
of multi-Terascale computers, the promise of full QCD calculations with light pion masses 
will finally be within our grasp. 

We note that subsequent to the completion of this work, a thorough analysis of the 
error correlations in the experimental measurements of polarized parton distributions was 



performed in ref |79|, with the result that the errors in the phenomenological values of 



moments of Am — Ad in table ITX| are of the order of ten percent. 
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APPENDIX A: CALCULATION OF SEQUENTIAL PROPAGATORS 

This appendix describes the explicit calculation of the sequential nucleon sources used 
in the calculations in this work. 

Denote the smeared field by the original unsmeared field by if) and the gauge invariant 
smearing function defined by Eq. [13] by F so that 

*o) = J d 3 xF^,(x ,x)ra>(x,t ). (Al) 

The forward propagator from a smeared source to a point sink, denoted 

S^,{^\x , t; x , t ) = t) f £(f„, t )), 

is the solution to the linear system 

dt' d 3 x' D w a S(x, t; x' , t)S£&(^\#, t'- x , t ) = 5(t - t ) F^(x, x ), (A2) 



where Dw denotes the Wilson Dirac operator and the smearing function is the source term. 
The propagator from a smeared source to a smeared sink is obtained by an additional 
smearing, 

S&{m\S, t; x , t ) d ^ (V a a (x, t) to)) 



d 3 x' x') S a a „ a a ,(mx\ t; xo, h). (A3) 

The three point function specifying u quark operator matrix elements, using the proton 



current, Eq. [12], with smeared quark fields, may be written: 

(J a OJ a ')( u \t h t s ,t ,Xi) = J d 3 x Q J d 3 x f e i ^e abc e a ' b ' c 'T (37 f (3 , Y x (A4) 
(U a a (x f , t f )U b (x f , tf)D%z f , t f ) u d v {x , t )OfX,(x> , U£(x h U)U b ;(x h U)D$fr, U)) 

where £j, tf and t Q are time coordinates of the source, sink and operator insertion, u and U 
denote point and smeared fields for the w-quark, d and D denote point and smeared fields 
for the (i-quark, and V = C75. The three point function for a d quark operator has the same 
form with the operator insertion d d ,(x" ,t )O dd ,,d d ,,(x l ,t' ). 

After performing the connected diagram contractions, the u and d operator matrix ele- 
ments may be written: 



(J a O J a ')^ d \t u t f ,t ) = J d 3 x J d 3 x f e^fM { ^ )[aa ' ]aa \x f ,t f ,x i ,t i ) (A5) 

x s«i(^\x f , t f] x , Q o d ss% (mz'o, c u) 

where 
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M { ^ [aa ' ]aa '(x f ,t f ,U) a ^ e abc e a ' b ' c ' S% 



def nbc^a'V 



X 



M^} aa>a \x fl t fl U)^ e abc e a ' b '' 



obb qcc , obb qcc 



(A6) 



(A7) 



We now define the backward propagator 

B [ ^'(x , t ,t f , a?,, U) d ^ J d 3 x f e-^f SX»(mto, t a ; x f , 2i,t f ) 7^ M*^' (x f , t f , U). 

(A8) 

Since it propagates from the sink to the operator, we use the relation / y 5 S(x, y)^ 5 = 
S^y, x) to obtain 



j d 3 Xf &** M [ ;f aa ' (x f , t f , x h u) siim\x f , t f - x 0l t c 



luu" nl (x ,t ,tf,Xi,ti) 



(A9) 



so that (|A5|) becomes 



(J a OJ a ')(t t ,t f ,t ) = J d 3 x o7 L"^v^ (A10) 

As in Eq. (|A2|) , the backward propagator B [aa ' ] defined by (|A8|) may be calculated by 
solving the linear system: 



dt d X Dw fj,"v"{p ) ^; %oito)$u» J {x >t ,tf) 

= 5(t-t f ) e~ 1 ^ f d 3 x' F a f (x, x)rf M*J;T' ]a ' >a ' (f, t f ). 



(All) 



Finally, in the same notation, the momentum projected two-point function with smeared 
sources is : 

(J a J a ')(t i ,t f ,x i ) = J d'xfe^fe^e^'T^T^ x (A12) 
([/«(£,, t f )U b (x f , tf)D;(x f , t f ) U£(x h U)V% 04 U)D c ; {x h U)) 



J d'xfe^fe^e^'V^f^S^, 



qaa' obb' cafe' qba' 

D aa' D l3f3 l D af3' D f3a' 
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TABLES 



I ~\ J J 

Data set 


QCD 


t 3 . , r 

x L t 


(3 


K>sea 




Approx. 
# configs. 


Irajectory 
separation 


SESAM 


J? 11 

full 


16^ x 32 


5.6 


0.1560 


0.1560 


200 


25 










0.1565 


0.1565 


200 


25 










0.1570 


0.1570 


200 


25 










0.1575 


0.1575 


200 


25 


SCRI 


full 


16 3 x 32 


5.5 


0.1596 


0.1596 


100 


20 










0.1600 


0.1600 


100 


20 










0.1604 


0.1604 


100 


18 


MIT 


quenched 


16 3 x 32 


6.0 


m q = oo 


0.1530 


200 














0.1540 


200 














0.1550 


200 




SESAM-cooled 


cooled 


16 3 x 32 


5.6 


0.1560 


0.1235 


100 












0.1570 


0.1246 


100 





TABLE I. Parameters specifying the full QCD, quenched, and cooled gluon configurations used 
in calculations of moments of quark distributions. For full QCD calculations, the number of hybird 
Monte Carlo trajectories between measurements is given in the last column. 
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observable H(4) mixing P lattice operator 



<4 a) 


6 3 + 


no 


1 


<4 6) 


3f 


no 







«r 


yes 


1 






* 

no 


1 


(1)a ? 

<*>2 


4 4 + 
63 
63 


no 

no 

no 

no 
no 




1 



1 








no 


1 


dl 




no** 





^2 


«r 


no** 


1 



4— 

<Z7{i£>4}<? 

qjW4q - ziqjiDiq + T1iD2q + TizDsq) 

qj{iDiD 4 }q ~ 2^(7{2-D2-D 4 } + 7{3£>3£>4})<7 

^{iDiDW^q + qi{2D2D3Ds}q - (3 <-> 4) 

<77 5 73<? 

<77 5 7{i-£>3}<? 

<?7 5 7{3-C>4}<7 

g7 5 7{i-D3-D4}9 
g7 5 cr 34 ^ 

97 5 ^3{4£>i}9 

97 5 7[3-D4]^ 

q7 5 7[iD { 3]D 4} q 



TABLE II. Operators used to measure moments of quark distributions. Different lattice opera- 
tors corresponding to the same continuum operator are denoted by superscripts a and b. Subscripts 
of irreducible representations of H(4) distinguish different representations of the same dimension- 
ality and superscripts denote charge conjugation C. In the operator mixing column, no* indicates 
a case in which mixing generically could exist but vanishes perturbatively for Wilson or over- 
lap fermions and no** indicates perturbative mixing with lower dimension operators for Wilson 
fermions but no mixing for overlap fermions. The entry in column P denotes the number of spatial 
components of the nucleon momentum, P, that must be chosen non-zero. Operators requiring one 
non-zero component have been written for P in the 1- direction and S in the 3-direction. 



observable 


7 


gLATT 




Z(/3 = 6.0) 


Z(J3 = 5.6) 




8/3 


-3.16486 


-40/9 


0.9892 


0.9884 




8/3 


-1.88259 


-40/9 


0.9784 


0.9768 




25/6 


-19.57184 


-67/9 


1.1024 


1.1097 


(* 3 ) q 


157/30 


-35.35192 


-2216/225 


1.2153 


1.2307 


(1)a 9 





15.79628 





0.8666 


0.8571 




8/3 


-4.09933 


-40/9 


0.9971 


0.9969 


<< 


8/3 


-4.09933 


-40/9 


0.9971 


0.9969 




25/6 


-19.56159 


-67/9 


1.1023 


1.1096 




1 


16.01808 


-1 


0.8563 


0.8461 


( X )Sq 


3 


-4.47754 


-5 


0.9956 


0.9953 


di 





0.36500 





0.9969 


0.9967 


d 2 


7/6 


-15.67745 


-35/18 


1.1159 


1.1242 



TABLE III. Renormalization constants 
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P = 6/sg 




ajv (fm) 


a p (fm) 


Ref. 


LANL 


5.6 


0.1585 


0.086(14) 


0.076(9) 


@ 




5.5 


0.16145 


0.116(7) 


0.104(5) 







5.4 


0.16450 


0.164(13) 


0.136(9) 





HEMCGC 


5.3 


0.16794 


0.124(4) 


0.108(2) 


m 


SCRI 


5.5 


0.16116 


0.109(4) 


0.090(3) 


m 


SESAM 


5.6 


0.1585 


0.091(16) 


0.076(10) 





TABLE IV. Parameters and lattice spacing for published dynamical Wilson fermion calcula- 
tions. This table summarizes calculations on 16 3 x 32 lattices and the resulting determinations of 



the lattice spacing from the masses of the nucleon, a at, 


and from the p meson 


a p . 




k c =0.1571 


0.1550 


0.1540 


0.1530 


<*>?' 


0.454(29) 
0.203(14) 


0.458(17) 
0.209(8) 


0.464(10) 


0.465(7) 




fl 30814 2M 
0.159(57) 


D 375(64"! 
0.179(29) 


0.402(47) 
0.178(21) 


0.437(41) 
0.192(19) 


\X ) u 

(x 2 )d 


0.0286(324) 


U.lO0(6Z ) 

0.0410(167) 


0.129(23) 
0.0464(118) 


0.142(20) 
0.0527(108) 


(X 3 }u 

(x 3 )d 


0.0369(357) 
-0.00893(1842) 


0.0441(192) 
0.000957(10126) 


0.0504(126) 
0.00917(638) 


0.0532(107) 
0.0132(53) 


(1)a« 
<l>Ad 


0.888(80) 
-0.241(58) 


0.915(46) 
-0.244(38) 


0.926(28) 
-0.250(18) 


0.939(20) 
-0.251(12) 




0.215(25) 
-0.0535(164) 


0.228(14) 
-0.0564(97) 


0.238(9) 
-0.0561(56) 


0.243(7) 
-0.0578(39) 


( x )a« 


0.141(123) 
-0.00144(7536) 


0.199(65) 
-0.0200(407) 


0.220(45) 
-0.0338(283) 


0.251(37) 
-0.0415(215) 


(a^Ad 


0.0269(428) 
0.00274(2530) 


0.0514(236) 
-0.00555(1431) 


0.0639(152) 
-0.0107(90) 


0.0758(121) 
-0.0145(66) 




1.014(82) 
-0.199(46) 


1.035(49) 
-0.222(27) 


1.045(28) 
-0.241(16) 


1.055(20) 
-0.251(11) 


(x)su 
{x)sd 


0.337(209) 
-0.0694(689) 


0.422(111) 
-0.0723(371) 


0.474(79) 
-0.0638(253) 


0.512(63) 
-0.0682(201) 


d\ 
df 


-1.354(75) 
0.282(42) 


-1.079(50) 
0.230(28) 


-0.914(24) 
0.208(13) 


-0.779(14) 
0.183(8) 


4 


-0.233(86) 
0.0396(311) 


-0.189(52) 
0.0313(193) 


-0.172(29) 
0.0230(104) 


-0.150(19) 
0.0193(65) 



TABLE V. Moments of quark distributions in quenched QCD calculated at three values of k 



and extrapolated linearly to k c . 
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k c =0.1585 


(0.1575) 


0.1570 


0.1565 


0.1560 




0.459(29) 


0.503(14) 


0.470(13) 


0.449(9) 


0.46 1(5) 


(x)i b) 

\ la 


0.190(17) 


0.222(9) 


0.207(8) 


0.207(4) 


0.214(4 ) 




462(1 74 s ! 


41 2(Q3l 


41 1 (691 


4SQ(70l 


D 394(451 




0.178(83) 


0.168(37) 


0.161(33) 


0.216(37) 


0.17 0(20) 


\X )u 

(x 2 ) r , 


U.l ( D^Do ) 

0.0314(303) 


n 1 qi (qa\ 

U.lol^o4 J 

0.0328(151) 


n 1 Q/i (OA \ 
0.0414(109) 


u.ioi^z4j 
0.0503(102) 


n 1 1 n/1 s "\ 
0.0496(101) 


(x 3 )i, 


0.0685(392) 


0.0443(168) 


0.0521(152) 


0.0594(121) 


0.0466(113) 


(x 3 ) d 


-0.00989(1529) 


0.0232(98) 


0.00789(549) 


0.0259( 63) 


0.0225(50) 


<1)a« 
(l)Ad 


0.860(69) 
-0.171(43) 


0.741(37) 
-0.214(25) 


0.880(28) 
-0.214(18) 


0.975(24) 
-0.248(15) 


0.936 (16) 
-0.254(9) 


<*>£! 


0.242(22) 


0.241(15) 


0.237(9) 


0.237(7) 


0.235(6) 




-0.0290(129) 


-0.0484(63) 


-0.0460(52) 


- 0.0605(38) 


-0.0621(34) 


(zA a) 


254(1 1 1 1 


205(831 


196(431 


21 7(41 1 


171 (311 

\J • J. f 11 ul J 




-0.0546(863) 


-0.0611(418) 


-0.0849(377) 


-0.0473(246) 


-0.0745(194) 


\X )Au 
{x 2 )Ad 


n ~\ 1R( AO\ 

0.00142(2515) 


u.uooy^oou j 
-0.0179(169) 


U.UD 1 O^IDI J 

-0.0149(101) 


u.uyzu^ioD J 
-0.0157(95) 


n HAR'il 1 1 A\ 
U.U4oo^ll4 J 

-0.0239(61) 


/1 \ 

<l>w 


U.yoo^oy ) 
-0.202(36) 


U.yly\4U J 
-0.239(27) 


-0.238(15) 


1 n^o/ 1 q\ 
l.Uoz^lo ) 

-0.225(14) 


1 PiTK /1 Q^l 
l.UfO (±o) 

-0.250(7) 




0.477(196) 
-0.144(68) 


0.424(109) 
-0.0828(385) 


0.418(79) 
-0.115(28) 


0.465(75) 
-0.0565(24 8) 


.405(48) 
-0.0739(167) 


d\ 
df 


-1.318(55) 
0.278(35) 


-1.036(44) 
0.269(30) 


-1.032(26) 
0.239(15) 


-0.957(15) 
0.200(13) 


-0. 854(11) 
0.197(6) 


4 


-0.228(81) 
0.0765(310) 


-0.191(50) 
0.0392(181) 


-0.179(34) 
0.0462(134) 


-0.164(29) 
0.0172(112) 


-0. 147(17) 
0.0154(59) 



TABLE VI. Moments of quark distributions in full QCD at ft = 5.6 using SESAM configura- 
tions. Lattice measurements are shown at four values of re, but to avoid finite volume effects at 
the lightest quark mass (re = 0.1575), only values at the three lowest re's are extrapolated linearly 
to re c . 
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/e c =0.16116 


0.1604 


0.1600 


0.1596 




0.416(41) 


0.459(25) 


0.451(16) 


0.480(9) 


\ id 


0.181(20) 


0.197(10) 


0.205(11) 


0.213(7) 


(r) (a} 


989(1 391 


379(68"! 


490(641 


D 469(471 


\ x )d 


U.U4oU(DZo ) 


nil 0(9^ 


n 1 89(90^ 
u. loz^zy j 


n 900(90^1 

u.zuy i^zu ) 


\X ) u 


U.Uiyi(4zy J 


n n7^Q/9i n\ 
U.U I Do(ZlU J 


n 1 ^nf QrA 
U.loU^oU J 


U. 140(1 1 ) 


(X 2 ) r , 
/ a 


0.0290(311) 


0.0435(171) 


0.0513(134) 


0.0590(94) 




0.0341(365) 


0.0401(201) 


0.0500(152) 


0.0509(112) 




0.0232(228) 


0.0229(138) 


0.00746(827) 


0.00975(583) 


(1)Au 


0.635(94) 


0.713(52) 


0.838(40) 


0.852(27) 


<l)Ad 


-0.279(48) 


-0.263(25) 


-0.280(22) 


-0.259(16) 




0.180(34) 


0.204(20) 


0.209(13) 


0.223(9) 


(^) Ad 


-0.0597(168) 


-0.0597(94) 


-0.0597(70) 


-0.0597(48) 


WAu 


278(1231 


217(661 


255(571 


192(391 




u. i iy ^oi ) 


n n/i i q( 31 

U.U4iy (OlO J 


-U.U14 1 ^ODO J 


-u.u4oy (zui ) 


\£ /A« 


U.1od(oU J 


u.uyo < (zoo J 


n 117/ 
U.ll ( (Zo J 


u.uoyo^ioo j 


\£ /Ad 


U.UU011(oU14 J 


-u.uzi4^ioy j 


U.UZoU^lOZ J 


-u.uzzo^yo j 


/1 \ r 
\l/<5u 


U.D / 0^ 1 0) 


U.o0o(4Z J 


u.you^oi j 


l.UoU^ZZ J 


\ l )Sd 


n 91 


-U.ZZO^OU J 


n 977(9^ 


n 9^^(14^ 

-U.ZDO^l^l J 


5u 


0.251(138) 


0.376(71) 


0.421(76) 


0.501(46) 


{x)sd 


0.00912(7595) 


0.0204(412) 


-0.106(31) 


-0.0527(253) 




-1.261(89) 


-1.197(56) 


-1.068(32) 


-1.046(21) 


d{ 


0.374(67) 


0.302(41) 


0.301(27) 


0.256(13) 


1U 

«2 


-0.264(77) 


-0.236(42) 


-0.203(40) 


-0.196(21) 


<g 


0.0191(402) 


0.00394(2320) 


0.0501(172) 


0.0266(100) 



TABLE VII. Moments of quark distributions in full QCD at (3 = 5.5 using SCRI configurations 



extrapolated linearly to k c . 
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k c =0.1266 


0.1235 


0.1246 




0.565(61) 


0.581(8) 


0.575(21) 


\ Id 


0.238(29) 


0.265(5) 


0.255(10) 


(r) (a} 


D 31 4(1 8Q1 


fl 555(74) 


n dP,R(dQ) 


\ x )d 


u. i iy ^oy j 


U.Z4 / ^OO J 


n 9ni ( r >A\ 

U.ZU1^Z4 ) 


\X ) u 


n 1 /i£/7c"\ 
U.140( (oj 


U.loy^zy J 


U.lol^Zl J 


(X 2 ) r , 
/ a 


0.0316(334) 


0.0716(123) 


0.0572(91) 




0.0517(337) 


0.0761(130) 


0.0673(88) 




0.00911(1390) 


0.0309(55) 


0.0231(35) 


(1)Au 


0.585(82) 


0.818(17) 


0.734(28) 


<l)Ad 


-0.298(46) 


-0.210(8) 


-0.242(16) 




0.118(46) 


0.114(17) 


0.115(13) 


(^) Ad 


-0.0120(242) 


-0.0230(57) 


-0.0191(79) 


(xA a) 


n 97Q( 3f0 


D 316(8") 


f) 303(9) 


/~,\( a ) 


-U.U / Do^ loo ) 


-U.U00U^40 J 


-u.uoy i ^oi ) 


\£ /Au 


U.ZoO^ll ( J 


n 07/1 { Ar\\ 
U.Z / 4^4U J 


n 9£i ( i >.A\ 
U.Zol^o4 J 


\£ /Ad 


-U.UDOU^ ( oZ J 


-U.U0o4^ZUU J 


-U.UD f O^ZOl ) 


/I \ r 


U. ( Do^DO ) 


u.yonio j 


U.ooD^ZZ J 


/1 \ r ■ 


-u.zoy^y j 




D 9^zlYl 7"\ 
-U.ZO^t^l / J 




0.451(265) 


0.702(96) 


0.611(73) 


(^)<5d 


0.0153(936) 


-0.137(27) 


-0.0820(288) 




-0.101(7) 


0.155(2) 


0.0628(22) 




0.00915(595) 


-0.0398(12) 


-0.0221(20) 


1U 

«2 


-0.0462(220) 


0.0619(95) 


0.0229(51) 


<g 


0.00673(1288) 


-0.0164(34) 


-0.00808(409) 



TABLE VIII. Moments of quark distributions calculated with cooled configurations, as de- 
scribed in the text, to eliminate most of the gluon degrees of freedom except instantons. Lattice 
measurements at k = 0.1235 and 0.1246 are extrapolated linearly to k c . 
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Connected 


QCDSF 


QCDSF 


Wuppertal 


Quenched 


Full QCD 


P henomenology 


M. E. 




(o = 0) 






(3 pts) 


(q±q) 


{x)u 


0.452(26) 






0.454(29) 


0.459(29) 




(x)d 


0.189(12) 






0.203(14) 


0.190(17) 




(x)u-d 


0.263(17) 






0.251(18) 


0.269(23) 


0.154(3) 


(X 2 )u 


0.104(20) 






0.119(61) 


0.176(63) 




(x 2 )d 


0.037(10) 






0.029(32) 


0.031(30) 




{X 2 )u~d 


0.067(22) 






0.090(68) 


0.145(69) 


0.055(1) 


(X 3 )u 


0.022(11) 






0.037(36) 


0.069(39) 




(x 3 ) d 


-0.001(7) 






0.009(18) 


-0.010(15) 




(X 3 )u~d 


0.023(13) 






0.028(49) 


0.078(41) 


0.023(1) 


(l)Au 


0.830(70) 


0.889(29) 


0.816(20) 


0.888(80) 


0.860(69) 




(l)Ad 


-0.244(22) 


-0.236(27) 


-0.237(9) 


-0.241(58) 


-0.171(43) 




(1}au-A<2 


1.074(90) 


1.14(3) 


1.053(27) 


1.129(98) 


1.031(81) 


1.248(2) 


{x)au 


0.198(8) 






0.215(25) 


0.242(22) 




(x)Ad 


-0.048(3) 






-0.054(16) 


-0.029(13) 




(x) Au—Ad 


0.246(9) 






0.269(29) 


0.271(25) 


0.196(9) 


(x 2 )au 


0.04(2) 






0.027(60) 


0.116(42) 




(x 2 )aj 

/Ad 


—0 012(6) 






-0 003(25) 


001 (25) 




(s 2 )A«-Ad 


0.05(2) 






0.030(65) 


0.115(49) 


0.061(6) 




0.93(3) 


0.980(30) 




1.01(8) 


0.963(59) 






-0.20(2) 


-0.234(17) 




-0.20(5) 


-0.202(36) 




<% 


-0.206(18) 






-0.233(86) 


-0.228(81) 




4 


-0.035(6) 






0.040(31) 


0.077(31) 





TABLE IX. Comparison of linear extrapolations of full QCD and quenched results with other 
lattice calculations and phenomenology at 4 GeV 2 in the MS scheme. The first column shows 
quenched results by the QCDSF collaboration at [5 = 6.0 [10,52-55] and the second column shows 
extrapolation of several moments to the continuum limit [15]. The third column shows full QCD 
results calculated using a different method with the same SESAM configurations we have used 
[14]. The quenched and full QCD results calculated in this work are shown in the fourth and 
fifth columns. Flavor non-singlet moments (x n ) of q{x) + (— l) n+1 q(x), Aq(x) + (— l) n Aq(x), and 
5q(x) + (—l) n+1 5q(x) are tabulated in the final column. Phenomenological unpolarized distributions 
are calculated from refs. [|]-[3| and polarized distributions are calculated from refs. [Q|| with error 
estimates as described in the text. 
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FIGURES 




FIG. 1. Connected (upper row) and disconnected (lower row) diagrams contributing to hadron 
matrix elements. The left column shows typical contributions of quarks and the right column shows 
contributions of antiquarks. 
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WS" MS" 

FIG. 3. Bootstrap distributions of {x)q for ensembles of 25 (left plot) and 204 (right plot) 
SESAM configurations at k = 0.1575 
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5.20 5.30 5.40 5.50 5.60 5.70 5.80 

P 

FIG. 4. The scaled pion mass squared, (am^) 2 x (a^^/a) 2 , for the dynamical Wilson fermion 
calculations tabulated in Table [fv|. The dotted line shows the mass at which the pion Compton 
wavelength equals one-fourth of the spatial dimension. 
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FIG. 5. Overlap between a smeared source and the proton ground state as a function of the 
source rms radius. The overlap for zero smearing is 6 x 10~ 5 . 




FIG. 6. Plateaus for (x)q calculated with point sources separated by 14 (left) and 12 (right) 
time steps using 100 configurations at j3 = 5.6 and k = .1575. 
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(a) 



(b) 




(c) 



t-t. 



(d) 



t-t. 



FIG. 7. Plateaus obtained using optimally smeared sources with a = 3 and JV = 40 separated 
by 12 time steps on SESAM configurations with (3 = 5.6 and k = 0.1560. Measurements as a 
function of Euclidean time of the operators (x)\ b \ (1)a<j, { x )^Aqi an< ^ (l)tfg are shown in panels a, 
b, c, and d respectively. Circles and squares denote matrix elements for up and down quarks 
respectively, the error bars are smaller than the symbols and shown within them, and the solid 
lines denote fits within the plateau region. 
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(x) 



FIG. 8. Comparison of operator measurements using combinations of point and smeared sources 
and sinks. From left to right, the four clusters, each containing four error bars, correspond to 
point-point, point-smeared, smeared-point and smeared-smeared source-sink combinations respec- 
tively. Within each cluster, from left to right, the four error bars correspond to windows in which 
the number of lattice points omitted from the window at the source and sink are (3,3), (3,5), (5,3), 
and (5,5) respectively. All smeared sources and sinks have N = 20 and a = 3. From top to bottom, 
the panels show measurements of the operators (x)^, (l)Ag, and (l)sq- 
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FIG. 9. Comparison of operator measurements using smeared sources of different sizes. From 
left to right, the three clusters, each containing four error bars, correspond to N = 0, N = 20, and 
N = 100. The windows within each cluster and the operators are the same as in Fig ^. 
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(x) 



V>) 



<1)A, 




(1)6, 




FIG. 10. Comparison of operator measurements using Dirichlet boundary conditions (left clus- 
ter) and periodic boundary conditions (right cluster). The windows within each cluster and the 
operators for the top three panels are the same as in Fig ||, and the bottom panel corresponds to 
the operator (x)a<j- 
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FIG. 11. Discrepancy between the axial charge, (V)Aq, calculated in single precision (solid 
curve) and double precision (dotted curve) as a function of (— log 10 where r^ n is the conju- 
gate gradient stopping residue. From the top panel to the bottom panel, the source-sink separation 
increases from 9 to 14 lattice spacings. 
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<x> 



.(b) 



— -fr — — 



•Quenched (MIT) 
Quenched (QCDSF) 



i-- 



0.2 



0.4 



0.6 



0.8 



m. 



7T 



GeV 



FIG. 12. Linear extrapolation of quenched calculations of the momentum fraction (x)^. The 
solid squares denote the results of the present calculation, and for comparison, QCDSF results are 
shown by the open points. Here and in subsequent figures, the upper and lower curves correspond 
to up and down quarks respectively. 
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FIG. 14. Linear extrapolation of quenched calculations of the axial charge, (1)a<j> where the 



symbols are as in Fig. |12. 
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FIG. 15. Linear extrapolation of quenched calculations of the moment (x)^ , where the symbols 
are as in Fig. 
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FIG. 16. Comparison of linear extrapolations of quenched and full QCD calculations of the 
momentum fraction, (x) q . The open symbols denote quenched calculations at f3 = 6.0 and the 
closed symbols denote full QCD calculations at (3 = 5.6. The upper and lower curves correspond 
to up and down quarks respectively. 
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FIG. 17. Comparison of linear extrapolations of quenched and full QCD calculations of the 
axial charge, (l)Aq, where the symbols are as in Fig. 
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FIG. 18. Comparison of linear extrapolations of quenched and full QCD calculations of the 
first moment of the quark spin distribution, (x)^ q , where the symbols are as in Fig. 16. 
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FIG. 19. Comparison of linear extrapolations of quenched and full QCD calculations of the 



tensor charge, (l)gq, where the symbols are as in Fig. 16 
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FIG. 20. Comparison of the quark momentum fraction (x)q^ calculated in full QCD and using 
configurations cooled to eliminate essentially all contributions except those of instantons. Solid 
symbols connected by solid curves denote full QCD and open symbols connected by dashed curves 
denote results after cooling. The upper and lower curves correspond to up and down quarks 
respectively. 
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FIG. 21. Comparison of the axial charge (1)aq calculated in full QCD and using cooled config 
urations, as in Fig. BO. 
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FIG. 22. Comparison of the first moment of the quark spin distribution, (x)^ q , calculated in 
full QCD and using cooled configurations, as in Fig. [20|. 
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FIG. 23. Comparison of the tensor charge (l)s q calculated in full QCD and using cooled con- 
figurations, as in Fig. [2(1. 
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FIG. 24. Comparison of the twist-three operator d\ from Eq. || calculated in full QCD and 
using cooled configurations, as in Fig. ^ . 
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FIG. 25. Chiral extrapolation of the momentum fraction using Eq. 21. Full QCD and 

quenched data calculated in this present work are denoted by squares and circles respectively, and 
the phenomenological result is indicated by the star. The least-squares fit and jackknife error bars 
are denoted by the heavy solid line and surrounding light lines. 
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FIG. 26. Chiral extrapolation of the first three moments of the proton quark distribution, 
(x n ) u _d, using Eq. £l[ Full QCD and quenched data calculated in this present work are denoted 
by diamonds and squares respectively, QCDSF quenched data [ 10 1 are denoted by triangles, and 
the phenomenological results are indicated by stars. The least-squares fit and jackknife error bars 
are denoted by the heavy solid line and surrounding light lines. 
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